# Analysis for Political Behavior article, "Economic Segregation and Unequal Policy Responsiveness".
# All analyses conducted using R version 3.5.1.


library(foreign)
library(lme4)
library(tidyverse)
library(effects)


# Change the file path below to directory where dataset is located.
setwd("/your/file/path/")



### Read data. ###

cces <- read.dta("CCES_merged.dta", 
                 convert.dates = TRUE, convert.factors = FALSE,
                 convert.underscore = TRUE, warn.missing.labels = TRUE)


# Baseline models.
# Dyadic representation.
mod1 <- lmer(formula = repscore ~ inc4 + pid.same +
                age + black + hisp + 
                expos8020 + zip.black.nh + zip.hisp + rural + 
                zip.gini + (1|zipcode) + (1|stdist) + (1|year), 
              data = cces, REML=F)
summary(mod1, correlation=FALSE)

# Dyadic representation measure without foreign policy issues.
mod1b <- lmer(formula = repscore2 ~ inc4 + pid.same +
               age + black + hisp + 
               expos8020 + zip.black.nh + zip.hisp + rural + 
               zip.gini + (1|zipcode) + (1|stdist) + (1|year), 
             data = cces, REML=F)
summary(mod1b, correlation=FALSE)

# Dyadic representation measure without foreign policy and social issues.
mod1c <- lmer(formula = repscore3 ~ inc4 + pid.same +
                age + black + hisp + 
                expos8020 + zip.black.nh + zip.hisp + rural + 
                zip.gini + (1|zipcode) + (1|stdist) + (1|year), 
              data = cces, REML=F)
summary(mod1c, correlation=FALSE)


# Campaign contacts.
mod2 <- glmer(formula = campcont ~ inc4 + 
                age + black + hisp + 
                expos8020 + zip.black.nh + zip.hisp + zip.gini + 
                rural + (1|zipcode) + (1|stdist) + (1|year), 
              data = cces, family = binomial(link="logit"), 
              nAGQ=0, control=glmerControl(optimizer = "nloptwrap"))
summary(mod2, correlation=FALSE)


# Donate.
mod3 <- glmer(formula = polact.donate ~ inc4 + 
                age + black + hisp + 
                expos8020 + zip.black.nh + zip.hisp + 
                zip.gini + rural + 
                (1|zipcode) + (1|stdist) + (1|year), 
               data = cces, family = binomial(link="logit"), 
               nAGQ=0, control=glmerControl(optimizer = "nloptwrap"))
summary(mod3, correlation=FALSE)


# Conditional effects by MC party affiliation.
# Create MC party dummy.
cces$mc.R[cces$party.H == 2] <- 1 # Rs
cces$mc.R[cces$party.H == 1] <- 0 # Ds

mod4a <- lmer(formula = repscore ~ inc4 + pid.same +
               age + black + hisp + 
               expos8020 + zip.black.nh + zip.hisp + rural + 
               zip.gini + (1|zipcode) + (1|stdist) + (1|year), 
             data = subset(cces, mc.R==1), REML=F)
summary(mod4a, correlation=FALSE)

mod4b <- lmer(formula = repscore ~ inc4 + pid.same +
                age + black + hisp + 
                expos8020 + zip.black.nh + zip.hisp + rural + 
                zip.gini + (1|zipcode) + (1|stdist) + (1|year), 
              data = subset(cces, mc.R==0), REML=F)
summary(mod4b, correlation=FALSE)


# Conditional effects by individual income.
# Effects package doesn't like variable names that end in numbers.
cces$inc <- cces$inc4

mod5 <- lmer(formula = repscore ~ pid.same +
               age + black + hisp + 
               expos8020 + inc + expos8020:inc +
               zip.black.nh + zip.hisp + rural + 
               zip.gini + (1|zipcode) + (1|stdist) + (1|year), 
             data = cces, REML=F)
summary(mod5, correlation=FALSE)

# Get segregation percentiles.
quantile(cces$expos8020, c(0.05, 0.25, 0.50, .075, 0.95), na.rm=TRUE)

eff.mod5 <- data.frame(effect("expos8020:inc", mod5, xlevels=list(expos8020 = c(-2.5, 1.5), inc = c(1, 2, 3, 4))))



### Alternative model specifications. ###

# Model year as factor/dummy.
# Dyadic representation.
mod7a <- lmer(formula = repscore ~ inc4 + pid.same +
               age + black + hisp + 
               expos8020 + zip.black.nh + zip.hisp + rural + 
               zip.gini + (1|zipcode) + (1|stdist) + factor(year), 
             data = cces, REML=F)
summary(mod7a, correlation=FALSE)

# Campaign contacts.
mod7b <- glmer(formula = campcont ~ inc4 + 
                age + black + hisp + 
                expos8020 + zip.black.nh + zip.hisp + zip.gini + 
                rural + (1|zipcode) + (1|stdist) + factor(year), 
              data = cces, family = binomial(link="logit"), 
              nAGQ=0, control=glmerControl(optimizer = "nloptwrap"))
summary(mod7b, correlation=FALSE)

# Donate.
mod7c <- glmer(formula = polact.donate ~ inc4 + 
                age + black + hisp + 
                expos8020 + zip.black.nh + zip.hisp + 
                zip.gini + rural + 
                (1|zipcode) + (1|stdist) + factor(year), 
              data = cces, family = binomial(link="logit"), 
              nAGQ=0, control=glmerControl(optimizer = "nloptwrap"))
summary(mod7c, correlation=FALSE)



### Tables. ###

library(stargazer)

stargazer(mod1, mod2, mod3, type = "text",
          star.cutoffs = c(0.05, 0.01),
          out = "table1.html")

stargazer(mod1, mod1b, mod1c, type = "text",
          star.cutoffs = c(0.05, 0.01),
          out = "tableA2.html")

stargazer(mod7a, mod7b, mod7c, type = "text",
          star.cutoffs = c(0.05, 0.01),
          out = "tableA1.html")

stargazer(mod4a, mod4b, type = "text",
          star.cutoffs = c(0.05, 0.01),
          column.labels = c("R subsample", "D subsample"),
          out = "tableA3.html")

stargazer(mod5, type = "text",
          star.cutoffs = c(0.05, 0.01),
          out = "tableA4.html")

